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We revisit the scenario of small-mass-ratio (g) black-hole binaries; performing new, more accurate, 
simulations of mass ratios 10:1 and 100:1 for initially nonspinning black holes. We propose fitting 
functions for the trajectories of the two black holes as a function of time and mass ratio (in the range 
1/100 < q < 1/10) that combine aspects of post-Newtonian trajectories at smaller orbital frequencies 
and plunging geodesies at larger frequencies. We then use these trajectories to compute waveforms 
via black hole perturbation theory. Using the advanced LIGO noise curve, we see a match of ~ 99.5% 
for the leading {£, m) — (2, 2) mode between the numerical relativity and perturbative waveforms. 
Nonleading modes have similarly high matches. We thus prove the feasibility of efficiently generating 
a bank of gravitational waveforms in the intermediate-mass-ratio regime using only a sparse set of 
full numerical simulations. 
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I. INTRODUCTION 

Numerical relativity (NR) has come a long way since 
the breakthroughs of 2005 [IHS] that allowed, for the first 
time, long-term evolution of black hole binaries (BHBs) . 
Among NR's significant achievements are its contribu- 
tions towards the modeling of astrophysical gravitational 
wave sources that will be relevant for the first direct de- 
tection and parameter estimation by gravitational wave 
observatories |4j . NR has also made contributions to the 
modeling of astrophysical sources, notably with the dis- 
covery of very large recoil velocities [IHS] , and the appli- 
cation of the numerical techniques to combined systems 
of black holes and neutron stars P . More mathematical 
aspects of relativity have also recently been investigated, 
including the evolution of N-black holes [10] , the explo- 
ration of the no- hair theorem pH, [T^ , and cosmic [T3] 
and topological censorship [T3], as well as BHBs in di- 
mensions higher than four [T5] . 

Among the remaining challenges are the exploration of 
the extremes of the BHB parameter space. The current 
state of the art simulations can simulate BHBs with mass 
ratios as small as g = 1/100 [121 HZ] and highly spinning 
BHBs with intrinsic spins a = Sh/M^ up to (at least) 
0.97 [TB]. Currently these runs are very costly and it is 
hard to foresee the possibility of completely covering the 
parameter space densely enough for match filtering the 
data coming from advanced laser interferometric detec- 
tors by the time they become operational. It is therefore 
imperative to develop interpolation techniques that allow 
for astrophysical parameter estimations, at a reasonable 
level of accuracy, based on a sparse set of numerical sim- 
ulations. 

In this spirit, we designed a set of prototypical runs 
for initially non spinning BHBs with small mass ratios 
q in the range 0.1 < q < 0.01. Since we expect that, 
for small enough mass ratios, this BHB system will be 



describable by perturbation theory, we compare the full 
numerical waveforms with those produced by perturba- 
tive evolutions via the Regge- Wheeler [19] and Zerilli f20] 
equations, supplemented by linear corrections in the spin 
of the large black hole [21] ■ The key ingredients pertur- 
bation theory needs is the relative trajectory of the small 
black hole with respect to the larger one and the back- 
ground mass and spin. In our previous tests we used 
the full numerical tracks and proved the perturbative 
waveforms agree reasonably well with the full numeri- 
cal ones [21] . In this paper we develop a model with free 
fitting parameters for these trajectories based on post- 
Newtonian and geodesic input, and fit to full numerical 
tracks. For these fits, we use full numerical evolutions of 
q = 1/10 and q = 1/100 BHBs [16, 21j and perform new, 
more accurate simulations. We compare the waveforms 
for the modeled tracks with the full numerical waveforms 
to confirm matching agreement within a fraction of a per- 
cent. Hence this method paves the way for further gen- 
eralizations and simulations to provide an approximate, 
yet accurate, bank for waveforms for second generation 
gravitational wave detectors. 



This paper is organized as follows. In Sec. |TT| we re- 
view the numerical methodology to perform the small 
mass ratio runs. In Sec. IIIII we describe the runs and 
results. In Sec. |IV| we describe the track modeling on 
post-Newtonian expansions and fits to the full numerical 
results. In Sec.|V|we give the results of using those tracks 
to generate perturbative waveforms and compares them 
with those extracted from the full numerical evolutions. 
We discuss the consequences and future application of 
this techniques in Sec. VI We also include an appendix [A] 
to briefly discuss perturbative theory in numerical coor- 
dinates (l-|-log, trumpet coordinates). 
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II. NUMERICAL RELATIVITY TECHNIQUES 

To compute the numerical initial data, we use the 
puncture approach [22] along with the TwoPuNC- 
TURES [53] thorn. In this approach the 3-metric on the 
initial slice has the form 7af, = {tpBL + u)'^Sab, where 
TpBL is the Brill-Lindquist conformal factor, Sab is the 
Euclidean metric, and u is (at least) on the punc- 
tures. The Brill-Lindquist conformal factor is given by 
ipBL = 1 + X]"=i ''iDj where n is the total num- 

ber of 'punctures', is the mass parameter of puncture 
i (m^ is not the horizon mass associated with puncture i), 
and fi is the coordinate location of puncture i. We evolve 
these black-hole-binary data-sets using the LazEv [23] 
implementation of the moving puncture approach [S] |3] 
with the conformal function W = ^Jx = exp(— 2(/)) sug- 
gested by Ref. [55]. For the runs presented here, we use 
centered, eighth-order finite differencing in space |10| and 
a fourth-order Runge-Kutta time integrator. (Note that 
we do not upwind the advection terms.) 

Our code uses the Cactus/EinsteinToolkit [5SII57] 
infrastructure. We use the Carpet 28^ mesh refinement 
driver to provide a "moving boxes" style of mesh refine- 
ment. In this approach refined grids of fixed size are ar- 
ranged about the coordinate centers of both holes. The 
Carpet code then moves these fine grids about the com- 
putational domain by following the trajectories of the two 
black holes. 

We use AHFinderDirect [21] to locate apparent 
horizons. We measure the magnitude of the horizon spin 
using the Isolated Horizon algorithm detailed in Ref. [30] ■ 
Note that once we have the horizon spin, we can calculate 
the horizon mass via the Christodoulou formula 



™2, + 5V(4mf„), 



(1) 



where mirr = ^/^/(IGTr) and A is the surface area of the 
horizon. We measure radiated energy, linear momentum, 
and angular momentum, in terms of ^/'4, using the for- 
mulae provided in Refs. [3TJ[32. However, rather than 
using the full V'4i we decompose it into t and m modes 
and solve for the radiated linear momentum, dropping 
terms with £ > 5. The formulae in Refs. [3T| (35] are 
valid at r = cx). Typically, we would extract the radi- 
ated energy-momentum at finite radius and extrapolate 
to r = oo. However, for the smaller mass ratios examined 
here, noise in the waveform introduces spurious effects 
that make these extrapolations inaccurate. We there- 
fore use the average of these quantities extracted at radii 
r = 70, 80, 90, 100 and use the difference between these 
quantities at different radii as a measure of the error. 

We extrapolate the waveform to r — oo using the per- 
turbative formula 1211 



hm [r^l'"(r,i)] 
= r^^(r,<)- 



+ 0(i?Obs)> 



(£-l)(^ + 2) 



where robs is the approximate areal radius of the sphere 
^Obs = const [Add a factor (1/2 — Mjr) multiplying the 
square bracket to correct for a difference in normaliza- 
tion between the Psikadelia (numerical) and Kinnersley 
tetrads at large distances.] We have found that this for- 
mula gives reliable extrapolations for i?obs lOOAf. 

Recent Cauchy-Characteristic extraction (CCE) stud- 
ies [33 (see also Ref. [53]) showed that this perturba- 
tive extrapolation formula can be more accurate than 
a linear extrapolation of the waveforms at finite radius 
to infinite resolution (provided that the observer i? is in 
the far zone). Those studies compared the extrapolated 
waveforms with the gauge invariant waveform on ob- 
tained using a nonlinear characteristic evolution from a 
finite radius to J^^ along outgoing null slices. The au- 
thors of Ref. |33j measured the errors in r?/'4 for an equal- 
mass BHB simulation when extracting sX R — 50Af and 
R = lOOAf, the corresponding extrapolation to 00, as 
well as the error in r'(/'4 obtained by applying the per- 
turbative extrapolation formula ([2]) to the R — lOOAf 
waveform. The errors in the perturbative extrapolation 
where the smallest over the entire waveform (both in am- 
plitude and phase). 

We note that multi-patch [35] and pseudospectral [36] 
techniques allow extraction radii very far from the source, 
leading to very small extrapolation errors. 



A. Gauge 

We obtain accurate, convergent waveforms and horizon 
parameters by evolving this system in conjunction with 
a modified l-flog lapse and a modified Gamma-driver 
shift condition [SJ [37], and an initial lapse a{t = 0) — 
2/(1 + ip%]^)- The lapse and shift are evolved with 

idt-l3'd,)a = ^2aK, (3a) 
= (3/4)^-77(x^^)/3^ (3b) 

where different functional dependences for ri[x°',t) have 
been proposed in Refs. [Ml [35M5] . Here we use a modi- 
fication of the form proposed in Ref. [39] , 



77(a;",t) = i?o 



(4) 



(2) 



where we chose Rq — 1.31. The above gauge condition 
is inspired by, but differs from Ref. [39] between the BHs 
and in the outer zones when a 7^ 1 and b ^ 2. Once 
the conformal factor settles down to its asymptotic ip = 
C / y'r + 0(1) form near the puncture, rj will have the 
form 77 = (i?o/C'^)(l -I- b{r/C'^)°') near the puncture and 
rj — Rqv''''^ M / (aM)^ as r — > 00. In practice we used 
a = 2 and 6 = 2, which reduces 77 by a factor of 4 at 
infinity when compared to the original version of this 
gauge proposed by Ref. [39] . We note that if we set 6 = 1 
then 7] will have a 1/r falloff at r = c» as suggested by 
Ref. [41j. Our tests indicate that the choices (a = 2, 
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TABLE I: Initial data parameters for the numerical simu- 
lations. Note that the q — 1/15 simulations are older and 
used a CFL factor twice as large. Here and rrij are the 
two puncture mass parameters, the puncture were located at 
(xijO, 0) and (a;2,0, 0) with momentum ±{Pr, Pt,0) and zero 
spin. The measured horizon masses and total ADM mass are 
also provided. 



Par am 


g= 1/10 


q = 1/15 


q = 1/100 


p 

m{ 


0.085237276 


0.057566227 


0.0086894746 


ml 


0.907396855 


0.936224183 


0.9896192142 


Xl 


7.633129115 


6.806172805 


4.952562636 


X2 


-0.7531758055 


-0.4438775230 


-0.04743736368 


Pt 


0.0366988 


0.0290721 


0.00672262416584 


Pr 


-0.000168519 


-0.000160518 


-0.00001026521884 


mui 


0.09129 


0.06254 


0.99065 


mH2 


0.91255 


0.94044 


0.00990841 


Maum 


1.00004 


1.00005 


1.00000000 


b = 1) and (a = 1,6 


— 1) lead to 


more noise in the 



0.0914 
0.09135 

0.0913 
0.09125 



n 0.0912 
E 



0.09115 
0.0911 
0.09105 
0.091 




h„ (dt/h=0.5) 
h„/1 .2 (clVh=0.5) 
h„/1 .2' (dt/h=0.5) 
h„ (dt/h=0.25) 
h„/1 .2 (dVh=0.25) 
h„/1 .2' (dVh=0.25) 



waveform than (a = 2, 6 = 2). 
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FIG. 1: The horizon mass conservation using CFL factors of 
0.5 and 0.25 for the q = 1/10 simulations. A factor of 10 
better results are obtained using a CFL factor of 0.25. The 
effect on the trajectory is significant because the mass loss 
is dynamically important (it significantly delays the merger). 
The mass changes post merger (sharp spikes near t ~ 950M 
and t ~ 1080M) are not dynamically important. 



III. SIMULATIONS AND RESULTS 

The initial data parameters for the q = 1/10, q = 1/15, 
and q = 1/100 simulations are given in Tablc|l] Note that 
the q — 1/15 simulations are older and suffer from the 
mass loss error discussed below. 

We used a base (coarsest) resolution of /iq = 4M with 
11 levels of refinement for q = 1/10 simulations and 15 
levels of refinement for q = 1/100 for the low resolu- 
tion simulations. The outer boundaries were at 400M. 
The higher resolution simulations were all based on these 
grids, but with correspondingly more gridpoints per level. 
The grid structure was chosen by studying the behavior 
of the background potential for the propagation of per- 
turbations [16' . In the appendix |A] we show this poten- 
tial, both in the isotropic coordinates of the initial data 
and final 'trumpet' coordinates. 

We previously evolved a set of g = 1/10, q = 1/15 [H], 
and q = 1/100 [TB^ BHBs using the standard choice of 
Courant-Friedrichs-Lewy (CFL) factor dt/h ~ 0.5. Al- 
though we found the waveform at lower resolution appear 
to converge, an unphysical mass loss led to incorrect dy- 
namics at later times (see Fig. [T]). This, in turn, led to 
oscillations in the errors as a function of grid resolution 
h. We found that reducing the CFL factor significantly 
reduces these unphysical effects [57]. In Figs, [l] and [2] 
we plot the horizon mass as a function of time for three 
resolutions using both the old and new CFL factor for 
both the q = 1/10 and q = 1/100 simulations. Note the 
much better conservation of the mass and the correspond- 
ing reduction in the lifetime of binary. Post-Newtonian 
trajectory evolutions (see Fig. Isl) indicate that the mass 
losses as small as 1 part in 10 are dynamically impor- 
tant, and the error introduced by this mass loss is sig- 
nificantly reduced by the new integration. The effects of 



the time integrator on the numerical error, and in par- 
ticular the mass loss and constraint violation errors, will 
be the subject of an upcoming paper by Ponce, Lousto, 
and Zlochower. 

For the q = 1/10 nonspinning BHB, we performed sim- 
ulations with central resolutions oi h = M/256, h = 
A//307.2, Af/337.9, M/368.64, and M/404.48, with 11 
levels of refinement in each case (the coarsest grid reso- 
lutions where hg, hQ/1.2, ho/l-i^, Hq/LM, ho/l.bS, re- 
spectively). We note that the waveform showed oscilla- 
tions in error as a function of resolution. We are therefore 
using the most widely spaced (in resolution) runs for the 
convergence plot. In particular, we use the h — M/256, 
M/337.9, and M/404.48 simulations. 

In Fig. [4] we show the convergence of the amplitude 
and phase of i/ja for these three resolutions. Convergence 
appears to break just at the start of the final plunge. 
Figure [5] shows the highest resolution waveform and the 
Richardson extrapolation (assuming a second-order er- 
ror). The extrapolation error becomes large just near the 
plunge, at a frequency of w = 0.19/M. We also plot the 
phase deference between the highest resolution run and 
the Richardson extrapolated waveform. The vertical line 
shows the point when the frequency is a; = 0.2/A/. The 
phase error, up to t = 850M is within 0.04 radian, but 
increases exponentially to 1.2 radians when uj — 0.2/AI. 

Determining convergence for the q — 1/100 simula- 
tions proved challenging for two reasons. There was a 
small random jump in the mass of the smaller BH be- 
tween resolutions that had a small effect on the trajec- 
tory. From Fig. [6] we can see a difference of one part 
m 10'' in the small BH mass. The two medium reso- 
lution have lower masses, and therefore merge slightly 
later than expected. The net effect is to confuse the or- 
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FIG. 2: The horizon mass conservation using CFL factors of 
0.5 and 0.25 for the q — 1/100 simulations. Here \5m\ is the 
absolute value of the difference between the expected horizon 
mass and the measured horizon mass for the smaller BH. In 
the case of g = 1/100, as in the case of g = 1/10, the horizon 
mass decreases (hence \Sm\ increases) with time when using a 
CFL factor of 0.5. A factor of 20 better results are obtained 
using a CFL factor of 0.25. The effect on the trajectory is 
significant because the mass loss is dynamically important (it 
delays the merger). The mass changes post merger (sharp 
spikes near t = 115Af and t = 172M) are not dynamically 
important. 



der when the mergers happens (we expected the higher 
resolutions runs to merge sUghtly later than the lower 
resolutions runs) . We see from Fig. [t] that the second- 
lowest resolution merges out of order with the other runs. 
Similarly, the third lowest resolution appears to merge 
late, as well. Another source of confusion for our con- 
vergence study has to do with lower-order errors (not 
necessarily associated with the mass) becoming impor- 
tant at higher resolutions. This leads to oscillations in 
the error as a function of resolution that can cause the 
merger time to oscillate with resolution. Nevertheless, 
the actual deviations in the trajectories, for the higher 
resolutions, is small. The improvement over the origi- 
nal q = 1/100 simulation is substantial. The mass con- 
servation is a factor of 20 better. The old simulations 
had an unphysical extra orbit inside the ISCO due to 
significant mass loss. Convergence measurements of the 
waveform do not appear to suffer much from the above 
mentioned mass loss. However, due to low amplitude of 
the waveform compared to the grid noise (induced by the 
refinement boundaries), the differences in the waveforms 
at the various resolution were smaller than the grid-noise 
fluctuations. Figures [9] and [lO] show the waveform and 
phase for the three highest resolutions. Note that the 
gridnoise apparent in Fig. |9] is due to reflections of the 
initial data pulse at the refinement boundaries. For our 
configuration, the imaginary part of the {£ — 2,m — 2) 
components of this pulse is small compared to the real 
part. This leads to a significant reduction (about a factor 
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FIG. 3: The 3.5PN orbital separations for the q = 1/10 bi- 
nary (the eccentricity is due to the fact that we use identical 
parameters as the full numerical simulation, which are not 
quasi-circular FN parameters) for slightly different small BH 
masses. A mass loss of 1 part in 10'* is dynamically impor- 
tant, leading to a significantly delayed merger. Reducing the 
mass change by a factor of 10 significantly reduces this error. 



TABLE II: Final remnant intrinsic spin, radiated energy, and 
recoil velocity for the q = 1/10, q = 1/15, and q — 1/100 
configurations. For reference we also include the predicted 
recoil velocity based on [43] . 



Par am 



a 

Vkick 
Vprcd 



g = i/io 



g = i/i5 



q = 1/100 



0.261 ±0.002 0.189 ± 0.006 0.0332 ± 0.0001 
0.0044 ± 0.0001 0.0022 ±0.0001 (5.5 ±1.0) x 10"^ 
60±2kms"* 34±2kms-* 1.0 ±0.1 km s"^ 
62 km s"* 34 km s"* 1.1 km s"^ 



of 10) in the noise in the imaginary part when compared 
to the real part. Figure [8] shows the trajectories for the 
four highest resolutions. Note that good agreement be- 
tween the curves when using the new smaller CFL factor. 

For reference, we report the final remnant properties 
for the three configurations in Table [ll] 

The above gauge conditions ([3| have a drawback near 
merger. The conformal function W goes to zero at the 
two punctures and there is a local maximum between 
the two punctures. Due to this local maximum, 77 ^ at 
some point between the two punctures. As the two punc- 
tures approach each other, this zero in r/ get progressively 
closer to the smaller BH. This, in turn, causes the coor- 
dinate radius of the horizon to shrink near merger. As 
seen in Fig. 11 the horizon radius for the new q — 1/100 
simulations decreases at merger, leading to a loss in res- 
olution of the smaller BH. Another issue related to this 
gauge is that the relative effective size (i.e. r^f/m^f) of 
the smaller BH is less than that of the larger BH. So the 
required resolution (gridpsacing) for each hole does not 
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FIG. 4: Convergence of the {£ — 2,m = 2) mode of i/)4 for 
the q — 1/10 waveform. The waveform starts out fourth- 
order accurate (due to our fourth-order waveform extraction 
algorithm), but then reduces to second order during the late 
inspiral, and dropping to first-order during the ringdown. The 
vertical line indicates the time when the frequency is tj = 
0.2/ M . The central resolutions for the three cases were hi = 
M/256, /i2 = M/337.9, and /13 = M/404.48, respectively, and 
each configuration used 11 levels of refinement. 



scale exactly with the BH's mass. This lower effective 
relative size of the smaller BH may be partially respon- 
sible for the reduced conservation of the horizon mass 
of the smaller BH when compared to the larger on (see 
Fig. p. 



FIG. 5: (TOP) A comparison of the highest resolution run 
with the Richardson extrapolated waveform. At the plunge, 
the phase error in the waveform is sufficiently large that the 
extrapolation give erroneous results (note the unphysical os- 
cillation near t = 1050M). The extrapolation breaks down at 
Lo = 0.19/M. (BOTTOM) The phase difference between the 
Richardson extrapolated waveform and the highest resolution. 



Schwarzschild perturbations. 

In the SRWZ formalism, the wave functions for the 
even and odd parity perturbations are expressed as a 
combined wave function of the first and second pertur- 
bative orders. 



it,r) + ^fj, {t,r), 



IV. PERTURBATIVE TECHNIQUES 



(o) 



(0,1) 



(o,Z,2) 
Ira 



(5) 



In Ref. HJ, we extended the Regge-Wheeler-Zerilli 
(RWZ) formalism [121 HO] of black hole perturbation the- 
ory to include, perturbatively, a term linear in the spin 
of the larger black hole. This Spin-Regge-Wheeler-Zerilli 
(SRWZ) formalism has second-order perturbations with 
linear dependence on the spin in the wave equations for 



where the superscript 1 and 2 denote the perturbative 
order, VP^m denotes the even parity Zerilli function, and 
^f^^ and ^-1^^'^^ denote the Regge- Wheeler (C unnmg- 
ham et al.) and the Zerilli functions for the odd parity 
perturbations (see Ref. [IJ for the relation between the 
wave functions). These functions satisfy the following 
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FIG. 6: The horizon mass for the q = 1/100 configurations. piG. 8: The puncture trajectories for the q = 1/100 configu- 
Jumps in the horizon mass of order 1/10^ are apparent be- ration, 
tween resolutions, this has a small but noticeable dynamic 
effect. 
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FIG. 7: The puncture separation as a function of time for 
the q = 1/100 configurations. The second lowest resolution 
run appears to merge later than expected, consistent with its 
reduced mass (apparent in Fig. |6| 



FIG. 9: The imaginary part (which is much less noisy than the 
real part) of the {£ = 2,m = 2) mode of i/'4 for the q = 1/100 
configurations. The differences in waveform between resolu- 
tions are dominated by grid noise at high resolutions. At the 
peak, the differences in the waveforms between resolutions is 
about 1% of the amplitude of the waveform. The real part of 
the waveform is roughly a factor of 10 more noisy. 
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FIG. 10: The phase of the {£ = 2,m = 2) mode of tp4 for the 
q — 1/100 configurations. The differences in phases between 
resolutions are dominated by grid noise. 
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FIG. 12: The horizon mass of both BHs for the new q = 1/100 
simulations. The horizon mass of the larger BH is conserved 
to a higher degree than that of the smaller BH. This is the 
case even though the number of gridpoints inside the smaller 
BH is actually larger than the number of gridpoints in the 
larger BH. 



three equations. 
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FIG. 11: The horizon radius of both BHs for the new q = 
1/100 simulations. The horizon radius of the larger BH barely 
changes, while the radius of the smaller BH decreases sharply 
at merger. 
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it, 



(6) 



where r* — r + 2Mln[r/{2M) — 1] is a characteris- 
tic coordinate, a is the non-dimensional spin parame- 
ter, y^(<=™"/°<i'i) denotes the Zerilli and Regge- Wheeler 
potential, respectively, Sgm is the source term, and 
(rp(t), ippit)) are the particle separation and orbital phase 
as a function of time. The differential operator p^°^™^ 
arises from the coupling between the BH's spin and the 
first-order wave functions 



We note that S'l^T''"^ and 



5, 



(odd,Z,2) 
In 



include local and extended source terms arising 



from the second perturbative order associated with prod- 
ucts of terms linear in a and the first-order perturbation 
functions. See Appendix A in Ref. [2T] for more details. 

In order to evolve the perturbative equations ([6]), we 
need the particle trajectory (rp(i), 0p(t)) as a function 
of time. Here we use a post-Newtonian inspired model 
for the BH trajectories, which will approximate the nu- 
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merical trajectories from intermediate-mass-ratio BHB 
inspirals. In order to model the numerical trajectories 
for the late inspiral phase of a BHB merger using a rel- 
atively straightforward fitting function, we start with an 
adiabatic quasicircular evolution, i.e., ignoring the eccen- 
tricity of the orbit and the component of radial velocity 
not due to the radiation reaction, and then parametrize 
deviations from a 3.5PN-TaylorT4 adiabatic inspiral. We 
model the time derivative of the orbital frequency dVt/dt 
and the orbital separation i? as a function of the orbital 
frequency 17 based on an extension of standard post- 
Newtonian calculations, described below. 

A. The orbital frequency Q. evolution of the NR 
trajectories 

First, we focus on the dfl/dt evolution with respect to 
the orbital frequency fl. In post-Newtonian (PN) cal- 
culations for the quasicircular case, the evolution for Q 
is obtained from the energy loss dE/dt and the relation 
between energy and frequency E{VL) given by 

dn _ dE (dE\^'^ 

~dt ^ ~dt \dn) ' ^ ^ 



where dE/dt and dE/dil are obtained by appropriate PN 
expansions. In the TaylorTl waveform, we simply sub- 
stitute these quantities into Eq. ([7|. By expanding the 
TaylorTl in the PN series again, we obtain the Taylor T4 
waveform, which has been shown to give good agreement 
with numerical simulations |45j . Therefore, we develop a 
fitting function based on the TaylorT4 evolution. These 
Taylor series waveforms are summarized in Ref. |46j . 



We note here that for Schwarzschild geodesies, dfl/dt 
diverges at the innermost stable circular orbit (ISCO), 
^Sch = 6A/ in Schwarzschild coordinates. This is be- 
cause dE/dVl becomes zero at the ISCO frequency Afil — 
(1/6)3/2 ^ 0.0680413817. Although we need to develop 
techniques for the transition to plunge [47l [48] in the 
above situation, there is no divergence in the Taylor T4 
waveform because the Taylor series expansion always has 
finite dVl/dt . 



When fitting of the trajectories in the NR coordinates 
for various mass ratios, we use the following modified 
Taylor T4 evolution. 



dVl 
~dt 
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(8) 



where d^l/dt and fl are obtained from the NR trajecto- 
ries, and M and ij denote the total mass and symmetric 
mass ratio of the binary system. 



M — nil + 1^2 ■ 
mim2 



(9) 



respectively. Here, A, a, B and fi are the fitting param- 
eters, and the power in the parameters should be a > 7 
and /3 > 7 in order to be consistent with this (7/2)PN 
formula. Here f2o is the frequency at i?sch ~ 3Af for 
circular orbits, i.e. M^a = il/if/'^ - 0.19245009. 

In Table |III[ we summarize the parameters obtained 
from fitting for the three cases, q = 1/10, q = 1/15 (from 



Ref. [21]), and q = 1/100. 

The fitting curve for the q — 1/10 case, which is valid 
up to Mn = 0.175, is shown in Fig. ^ The fitting 
parameters, a and /3, for this case are consistent with the 
current PN formula, and we see that the fitting curve is 
close to the TaylorT4 evolution for small MVl. 

We show the fitting curve for the q— 1/15 case, which 
is vahd up to MQ. = 0.17, in Fig. [it] Again, the fitting 
parameters, a and /3, for this case are consistent with the 
PN formula. 



For the q — 1/100 case, which is shown in Fig. 15 



we 



first found by a direct fit of all parameters that a < 7, 
which is not consistent with the PN approach. One pos- 
sibility is that the trajectories are too short. If we were to 



TABLE III: The summary of the fitting parameters for the 
"PN+" (a modification of the Taylor T4 evolution) in Eq. 
The first line for the q = 1/100 case (*) gives the fitting 
parameters obtained by allowing all 4 constants to vary, the 
second line (f) for the q — 1/100 case gives the parameters 
from a fitting of A and B only by assuming a = 8 and 13 — 15, 
and the third line (J) gives the parameters extrapolated from 
the q = 1/10 and q = 1/15 results using Eq. (111. 



Mass ratio 


A 


a 


B 




g = l/10 


17.0500 


7.21975 


8.18920 


12.5197 


g = l/15 


26.0150 


7.54047 


8.65525 


13.6168 


g= 1/100* 


93.0650 


4.32071 


5.42457 


14.9711 


q = 1/1001" 


288.888 


8 


16.2033 


15 


q = 1/100* 


233.985 


9.45214 


11.5411 


21.0731 



FIG. 13: The fitting curve for the g = 1/10 case using Eq. g 
and orbital frequencies up to MQ = 0.175. The solid (black), 
dashed (red) and dotted (blue) curves show the NR, fitting 
and TaylorT4 PN evolutions, respectively. In the NR evolu- 
tion, the end point of the solid (black) curve corresponds to 
JJsch = 2M, and MQ = 0.175 is the NR orbital frequency 
around -Rsch = 3M, assuming the NR coordinate system can 
be approximated by a "trumpet" stationary 1 + log slice of 
the Schwarzschild spacetime. The inset shows the zoom-in of 
the initial part (around MJl = 0.05) 




extend them to lower frequencies, we would expect that 
the slope increases, leading to a larger a. Below, we try 
two alternate fitting methods in order to try to improve 
the fit. This first fit gives the best agreement with the 
numerical trajectory so this is the one we use in Sec. [V| 
to calculate the gravitational waveforms. However, the 
second fit reproduces the expected behavior at large sep- 
arations, so we expect that it gives a better waveform in 
that regime. 

Note that, as seen in Fig. [15] the end point of the 
NR trajectory, which corresponds to i?sch = 2A/, has 
dfl/dt < 0. This is a feature of Schwarzschild geodesies in 
Schwarzschild coordinates (e.g., a quasi-circular geodesic 
sfightly inside the ISCO has dn/dt = -0.0113). Hence 
it appears that in the q = 1/100 case the trajectories are 
much closer to geodesic motion than in the other cases. 

Due to the fact that the q = 1/100 simulations were 
very short, obtaining accurate fitting parameters for the 



FIG. 14: The fitting curve for the g = 1/15 case using Eq. Q 
up to orbital frequencies of Mfl — 0.17. The solid (black), 
dashed (red) and dotted (blue) curves show the NR, fitting 
and TaylorT4 PN evolutions, respectively. In the NR evo- 
lution, the end point of the solid (black) curve corresponds 
to -Rsch ~ 2M, and Mfl = 0.17 is the NR orbital frequency 
around -Rsch ~ 3M, assuming the NR coordinate system can 
be approximated by a "trumpet" stationary 1 + log slice of 
the Schwarzschild spacetime. 
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FIG. 15: The fitting curve for the g = 1/100 case using 
Eq. ([8| with orbital frequencies up to AdQ — 0.15. The solid 
(black) and dotted (blue) curves show the NR and TaylorT4 
PN evolutions, respectively. The Fitting (*, dashed (red)), 
(1, dot-dashed (green)) and (2, dot-dot-dashed (magenta)) 
curves show the fitting curve with the (*), (f ) and (J) cases in 
Table [ml In the NR evolution, Affi = 0.15 is the NR orbital 
frequency around -Rsch = 3M, assuming the NR coordinate 
system can be approximated by a "trumpet" stationary 1+log 
slice of the Schwarzschild spacetime. 
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trajectory is error prone. We therefore consider several 
different techniques. If we assume a simple fitting with 
respect to the mass ratio, 



{A, a, B, /?) = c^jf^ , 



(10) 



for the fitting parameters in the case of the q — 1/10 and 



10 



1/15 cases, we obtain 



A = 0.797021 ?7~i-22855 

a = 5.26843 7?-°i26379^ 

B = 5.48252 77-°i6°938^ 

/3 = 6.80971 77-"-244245_ 



(11) 



Then, inserting the symmetric mass ratio 77 — 10/121 for 
q = 1/100 gives fitting parameters for q — 1/100. We 
report these values on the third hne labeled q = 1/100 
in Table |llT] The curve obtained by these parameters is 
denoted by "Fitting (2)" in Fig. 15 Note that these pa- 



rameters actually produce a relatively poor fit for the 
trajectory. On the other hand, if we assume a — 8 and 
/3 = 15 and fit the remaining parameters using the NR 
trajectory, we obtain the "Fitting (1)" parameters (sec- 
ond hne labeled q = 1/100 in Table[in|. This fit smoothly 
interpolates between the PN trajectory at small frequen- 
cies and the start of the numerical trajectory at higher 
frequencies. Since both the standard fitting procedure. 



and "Fitting(l)" give reasonable fits, but different pre- 
dictions, in order to accurately determine the constant in 
Eqs. ( [IT] ), we need both more accurate and longer evolu- 
tions NR simulations of g = 1/100. 



B. The orbital radius R versus orbital frequency Q, 
in the NR coordinates 



In order to complete our description of the trajectories 
we need a fitting function that will allow us to express 
the separation i? as a function of Q,. We use a func- 
tion R{n) inspired by the ADM-TT PN approach gS], 
because the isotropic coordinates (used in the ADM-TT 
PN approach) and the radial isotropic "trumpet" station- 
ary 1 -|- log slices of the Schwarzschild spacctime are very 
similar [50l I51j . 

The fitting function based on the ADM-TT PN ap- 
proach is given by 



R 



M 



1 



(Aff])2/3 

1 1625 



-1 + ^r;) (Amf/3 



N4/3 



1 67 2 

v + T7^v^'-ov' + ^v'] {Mnf ) / (1 + a,{n/nor) + c , 



192 



81 



(12) 



TABLE IV: The fitting parameters for the "PN-I-" (a modi- 
fication of the ADM-TT PN calculation) in Eq. ([l2| for the 
relation between R and 57 in the NR coordinates. 



Mass ratio 


C 


ao 


ffli 


g=l/10 


0.216953 


0.513214 


4.68472 


g=l/15 


0.237427 


0.600321 


4.57899 


q = 1/100 


0.198137 


0.923360 


5.29681 



where R and are in the NR coordinates, and ao, ai and 
C are fitting parameters, and A/ilg = (1/3)'^/^. Note that 
ai should be greater than 2 to be consistent with the PN 
calculation. 

Initially, we expected that C, which is the difference 
between the NR radial coordinate and "trumpet" coordi- 
nate, would be very small because this term is inconsis- 
tent with the ADM-TT PN formula. However, we found 
that this term is non-trivial for all three mass-ratio cases. 

The results of the fits are shown in Figs. [T6]|T8| for 
q — 1/10, q = 1/15, and q = 1/100, respectively, and the 
fitting parameters are summarized in Table |IV| We see 
from Table [IV| that the values of ai for all cases are con- 
sistent with the current PN calculation (ai > 2), and find 
a non-negligible C ~ 0.2 for all three cases. Due to this 
C contribution, there is a finite difference between the or- 



bital radius evaluated in the ADM-TT PN approach and 
that by the fitting function in Eq. (12) even for smaller 
MVL. 

As a result of the above fitting analysis, it is necessary 
to consider a radial transformation to remove the offset 
C between the NR coordinates (i?NR — R for Eq. (12)) 
and the "trumpet" coordinates, 



-Rnr, Rhog ~ Rnr — C . 



(13) 



After removing this offset, we then transform to the stan- 
dard Schwarzschild coordinates, where the explicit time 
and radial coordinate transformations, 

(TLog, -RLog) — > (Tsch, ^Sch) , (14) 

are given in Rcf. [52 . Here, we assume Tnr = TLog- We 
note that after applying the C transformation discussed 
above, we find that the frequency when i?sch — 3Af is 
given by Mn = 0.161, 0.158 and 0.152 for the q ^ 1/10, 
q = 1/15 and q = 1/100 cases, respectively. 

Note that this procedure of fitting NR trajectories di- 
rectly in the NR coordinates and then transforming to 
Schwarzschild coordinates, appears simpler than the al- 
ternative of evolving the perturbations directly in the 
"trumpet" coordinate system. We nevertheless study 
this alternative in the Appendix |Aj where we provide the 
sourceless part of the perturbative equation and study its 
characteristic speeds in "trumpet" coordinates. 
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FIG. 16: The fitting of the orbital radius R vs. orbital fre- 
quency Slin the NR coordinates for the q = 1/10 case. The 
solid (black), dashed (red) and dotted (blue) curves show the 
NR, fitting and ADM-TT PN evolutions, respectively. The 
fitting by using Eq. ([12} is vahd up to Mil = 0.175. In the NR 
evolution, the end point of the solid (black) curve corresponds 
to Rsch = 2M, and MQ = 0.175 is the NR orbital frequency 
around -Rsch = 3M, assuming the NR coordinate system can 
be approximated by a "trumpet" stationary 1 + log slice of 
the Schwarzschild spacetime. 




FIG. 18: The fitting of the orbital radius R vs. orbital fre- 
quency Q in the NR coordinates for the q = 1/100 case. The 
solid (black), dashed (red) and dotted (blue) curves show the 
NR, fitting and ADM-TT PN evolutions, respectively. The 
fitting by using Eq. Q is vahd up to MQ = 0.15. In the NR 
evolution, the end point of the solid (black) curve corresponds 
to Rsch = 2M, and MQ, = 0.15 is the NR orbital frequency 
around -Rsch ~ 3M, assuming the NR coordinate system can 
be approximated by a "trumpet" stationary 1 + log slice of 
the Schwarzschild spacetime. 



NR(q=l/IOO) 

— - Filling 
■ ■ . PN 



FIG. 17: The fitting of the orbital radius R vs. orbital fre- 
quency Q in the NR coordinates for the q = 1/15 case. The 
solid (black), dashed (red) and dotted (blue) curves show the 
NR, fitting and ADM-TT PN evolutions, respectively. The 
fitting by using Eq. (p2| is vahd up to MQ = 0.17. In the NR 



evolution, the end point of the solid (black) curve corresponds 
to i?sch = 2M, and MQ = 0.17 is the NR orbital frequency 
around -Rsch = 3M, assuming the NR coordinate system can 
be approximated by a "trumpet" stationary 1 + log slice of 
the Schwarzschild spacetime. 




V. RESULTS 

We obtain an approximation to i7(TNR) by integrat- 
ing Eq. We then obtain i?(rNR) using Eq. ^ and 
the approximate ^(Tnr). The initial orbital phase is 
not fixed, and we set it here such that the matching 
between the perturbative and numerical waveforms are 
maximized. 



Note that part of the difference between the NR wave- 
form and the perturbative waveforms is due to eccentric- 
ity in the numerical trajectories. This affects the time 
step of the orbital evolution through 



dt 



dil. 



(15) 



This may help explain why in the q = 1/10 case, we 
observe a slower time evolution in the fitting trajectory 
than the NR evolution. 

In order to generate the final plunge trajectory and 
regularize the behavior of the source near the horizon in 
the SRWZ formalism [3T], we attach the model for the 
trajectories to a plunging geodesic in a continuous way. 
To do this, we choose a matching radius Rm ~ 3M in 
Schwarzschild coordinates (the light ring) and then de- 
termine the energy and angular momentum of the tra- 
jectory at that point. We then complete the trajectory 
from i? ~ Rm to i? ~ 2M using the geodesic with the 
same energy and angular momentum. 

When we set the matching radius to Rm = 3M, we 
find that the energy E and angular momentum L in the 
Schwarzschild coordinates are given by 



E 
L 



3"'. 

9 M 2 u'f 



(16) 



where u'^ = dx^ /dr is the four velocity, and we focus on 
equatorial orbits. To evaluate u* in the above equation, 
we use 



-1. 



(17) 
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where g^i, denotes the Schwarzschild metric, and at 
Rm = 3M this becomes 



- 3(i?sch(Tsch))' 



-9M2($sch(rsch))' 



-1/2 



(18) 



Here, i?sch = w'^/w* = di?sch/dJsch and $sch = u"^/"* 
d^Sch/dTsch are the velocity components obtained from 
the numerical data after the coordinate transformation 
from the "trumpet" coordinat es t o the Schwarzschild co- 
ordinates. We also use Eq. (17) to calculate in the 
source term of the SRWZ formalism. 

In the SRWZ formalism, the gravitational waveforms 
at a sufficiently distant location i?obs is given by 



2M 



We use the normalized waveform, i.e., the expression in 
the right hand side of the above equation to calculate the 
wave amplitude. 

On the other hand, to convert the NR ip4, which is ex- 
trapolated to infinity via Eq. ([2]), to the waveform h, we 
use the "pyGWAnalysis" code in EinsteinToolkit [27^ 
(see Ref. [53] for details). Then, we calculate the match 
between the NR and SRWZ waveforms in the advanced 
LIGO (Zero Det, High Power) noise curve [M]. In Ta- 
blc|Vj wc summarize the match for the three mass ratios. 
In the following subsections, we treat each mass ratio case 
separately. 



TABLE V: The match between the NR and SRWZ waveforms 
for various total masses using the advanced LIGO noise curve. 
For the smaller masses, the entire waveform is in the advanced 
LIGO band, while for the larger masses, only the final ring- 
down part of the waveform is in band. For these larger masses, 
we integrate over frequencies / > lOHz. Here we only used 
the {£ — 2, m = 2) mode to calculate the match. 



Total mass (Mq) 


g= 1/10 


g= 1/15 


q = 1/100 


100 


0.983019 


0.982900 


0.995162 


200 


0.992766 


0.992704 


0.995216 


300 


0.996219 


0.996514 


0.995486 


400 


0.996781 


0.997734 


0.995513 


500 


0.996844 


0.997973 


0.995474 



A. The g = 1/100 case 

For the q = 1/100 case, the fitting formulas in Eqs. ([s]) 
and ( 12 ) are valid for orbital frequencies up to Mfi = 0.15 
in the "trumpet" coordinates. Here, we extend them 



beyond Mf2 = 0.15 by attaching these trajectories to a 
plunging geodesic at Rm (as explained above). 

Ideally, we would like to use as much of the fitting tra- 
jectory as possible. However, when attempting to use 
the fitting trajectory to calculate the 4-velocity at small 
radii, we can not evaluate u' for i?sch < 2.12M from 



Eq. ( 17) because the approximation that the background 
metric is Schwarzschild in "trumpet" coordinates breaks 
down and the vector (1, i?sch, 0, $sch) becomes spacelike 
in the background metric (but not in the numerical met- 
ric). This breakdown can be thought to arise from the 
interaction of the singular part of the conformal factor 
due to the smaller BH (when it is close to horizon of the 
larger one) with the singular part of the conformal fac- 
tor due to the larger BH. This, in turn, changes and 
therefore the gauge, when the two BHs approach each 
other. Also, since we see some delay in the phase evolu- 
tion for a small matching radius near the horizon, we set 
Rm ^ 3M for the start of the geodesic approximation. 
At this matching radius, E and L for the Schwarzschild 
geodesic are evaluated from the orbital separation and 



three velocity using Eq. (16). We find 



E = 0.93942436, 
L/AI = 3.46184226. 



(20) 



In the SRWZ formalism, we set the non-dimensional spin 
parameter a — 0.033 (see Table [lT|) . Although this value 
is obtained in the NR simulation, we can also obtain it 
from an empirical formula 55]. 

We next calculate the match in the advanced LIGO 
noise curve. For the integration in the frequency domain, 
we consider gravitational wave frequencies Mfl22 > 0.15. 
For a total mass M = 484M0, where we have the lower 
frequency cut off (corresponding to the start of the nu- 
merical waveform) at /low = lO.OlHz, the match between 
the NR and SRWZ waveforms is 



M22 = 0.995477, 



(21) 



for the {£ — 2, m = 2) mode. The amplitude and phase 
evolutions for the NR and SRWZ waveforms are shown 
in Fig. 



19 



To obtain this figure, we translated the wave- 
forms so that the maximum in the amplitude occurs at 
t ~ and used the freedom in the phase to set it to zero 
at t = -75M. 

The {£ — 2, m — 1) mode of h is the leading contri- 
bution from the odd parity in the RWZ calculation. The 
amplitude and phase evolutions for this mode are shown 
in Fig. 20 After adding a time translation so that the 
maximum in the amplitude occurs at t = and a phase 
translation so that the phase 4> — is zero a.t t — ~75Af , 
we still see a small difference between the phases. For a 
total mass of M = 242M0 (this mass is chosen such that 
the lowest frequency is just in the advanced LIGO band 
= lO.OlHz), the match between the NR and SRWZ 
waveforms for this mode is given by 



M2 



0.957966. 



(22) 
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FIG. 19: The amplitude (TOP) and phase (BOTTOM) of 
the {£ = 2,m = 2) mode of h for the q = 1/100 case. The 
waveforms were translated in time such that the maximum in 
the amplitude occurs at t — and the phases were adjusted 
by a constant offset such that (p = at t = —75M. The 
solid (black) and dashed (red) curves show the NR and SRWZ 
waveforms, respectively. 
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FIG. 20: The amplitude (TOP) and phase (BOTTOM) of 
the {£ = 2, m = 1) mode of h for the q — 1/100 case. The 
waveforms were translated in time such that the maximum in 
the amplitude occurs at t — and the phases were adjusted 
by a constant offset such that tf) = a.t t — —75M. The 
solid (black) and dashed (red) curves show the NR and SRWZ 
waveforms, respectively. 
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For the {£ ^ 3, m 



phase evolutions are shown in Fig. 21 
(mass chosen so that /low 
mode is 



3) mode of h, the amplitude and 
For M = 726M0 
the match for this 



lO.OlHi) 



M33 = 0.993895. 



B. The (7 = 1/10 case 



(23) 



For the q = 1/10 case, the fitting formulae is valid up 
to larger orbital frequencies than for the q = 1/100 case. 
Nevertheless, u* can not be evaluated for i?sch < 2.26Af 
from Eq. (17). Therefore, we need to start the geodesic 



approximation at radii larger than this, and again we set 
the matching radius i?M ~ 3M. We find that the energy 
and angular momentum for the Schwarzschild geodesic 
are given by 



E = 1.02546489, 
L/M = 3.95926146, 



(24) 



at the matching radius. Not that the energy is £' > 1 
because we simply project the orbital quantities on the 
Schwarzschild spacetime by using Eqs. ( 16 ) and (171 (this 
is a strong indication that the motion of the q = 1/10 
binary is non-geodesic even close to the larger BH). In 
the SRWZ formalism, we set the non-dimensional spin 
parameter to a = 0.26, consistent with the measured 
spin from the numerical simulation (see Table |ll| . 

The amplitude and phase evolutions of the {£ = 2, m = 
2) mode of h are shown in Fig. 22 While the phases of 



the NR and SRWZ waveforms are in good agreement, 
there are relatively large differences in the amplitudes 
both for the inspiral and merger phases. We will discuss 
the amplitude differences in more detail in Sec. [VD| 

We consider the match calculation for MVl.22 > 0.075 
and total mass M = 2A2Mq, i.e., the lower frequency cut 
off at /low — lO.OlHz. In this case the match between 
the NR and SRWZ waveforms is 



M-2 



0.994669 . 



(25) 



for the (€ = 2, m = 2) mode in the advanced LIGO noise 
curve. 
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FIG. 21: The amplitude (TOP) and phase (BOTTOM) of 
the {£ = 3, m — 3) mode of h for the q — 1/100 case. The 
waveforms were translated in time such that the maximum in 
the amplitude occurs at t — and the phases were adjusted 
by a constant offset such that (p = at t = —75M. The 
solid (black) and dashed (red) curves show the NR and SRWZ 
waveforms, respectively. 




FIG. 22: The amplitude (TOP) and phase (BOTTOM) of 
the (^ = 2, m = 2) mode of h for the q — 1/10 case. The 
waveforms were translated in time such that the maximum in 
the amplitude occurs sd, t — and the phases were adjusted 
by a constant offset such that = at i = — 830M. The 
solid (black) and dashed (red) curves show the NR and SRWZ 
waveforms, respectively. 




The amplitude and phase evolutions for the {£ 
1) mode of h are shown in Fig. 23 , For M 



m 



I2IM0 (/low = lO.OlHz), the match for this mode is 



M2 



0.989382 . 



(26) 



The amplitude and phase evolutions for the {£ = 
3, TO = 3) mode of h are shown in Fig. 24 For M = 
363M0 (/low — lO.OlHz), the match for this mode is 

M33 = 0.987414. (27) 

We have a similar mismatch in the (£ = 2, m = 1) and 
(£ = 3, TO = 3) modes. 

We note that although the = 2, to = 1) and {£ — 
3, TO = 3) modes have a slightly large mismatch (1 — 
Ai), these modes are the sub-dominant contributions to 
the gravitational wave. The maximum amplitudes are ^ 
0.12, 0.035 and 0.04 for the (£ = 2, to = 2), {£ = 2, m = 
1) and (^ = 3, TO = 3) modes, respectively. Therefore, 
we find that from a rough estimation the contribution 
of the mismatch for each mode are almost same in the 
gravitational wave. 



Also, as discussed in Ref. [ST, unlike the case of q = 
1/100, here the spin corrections in the SRWZ formalism 
are important for obtaining the correct ringdown phase. 
We see that the correct ringdown frequency was obtained 
by comparing the slope of the phase evolution in Figs. [22] 
[21] and [Ml 



C. The q = 1/15 case 

For the q = 1/15 case, the fitting formula is valid up 
to larger orbital frequencies than for the q = 1/100 case. 
Here too, we can not calculate w* for i?sch < 2.21Af from 
Eq. (17), and again we use the geodesic approximation 

We find that the 



and a matching radius Rm ^ 3M. 
energy and angular momentum are given by 

E = 0.985372666, 
L/M = 3.74699937. 



(28) 



In the SRWZ formalism, we set the non-dimensional spin 
parameter a = 0.19, consistent with the numerical results 
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FIG. 23: The amplitude (TOP) and phase (BOTTOM) of 
the {£ = 2, m = 1) mode of h for the q = 1/10 case. The 
waveforms were translated in time such that the maximum in 
the amplitude occurs at t — and the phases were adjusted 
by a constant offset such that (j> = at t — — 830M. The 
solid (black) and dashed (red) curves show the NR and SRWZ 
waveforms, respectively. 




FIG. 24: The amplitude (TOP) and phase (BOTTOM) of 
the (^ = 3, m = 3) mode of h for the q — 1/10 case. The 
waveforms were translated in time such that the maximum in 
the amplitude occurs sd, t — and the phases were adjusted 
by a constant offset such that = at i = — 830M. The 
solid (black) and dashed (red) curves show the NR and SRWZ 
waveforms, respectively. 




(see Table |n]). 

The amplitude and phase evolutions of the 
2) mode of h are shown in Fig. 25 
Mn22 > 0.09 with a total mass of M 
10.02Hz) is 

M22 = 0.996039, 



= 2, rn = 
The match for 
29OA/0 (/low = 

(29) 



using the advanced LIGO noise curve. 
The amplitude and phase evolutions 



for 
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the {£ 
For M 



2, TO = 1) mode of h are shown in Fig, 
145M0 (/low = 10.02Hz), the match for this" mode is 

M21 ^ 0.989005. (30) 

The amplitude and phase evolutions for the {£ = 



3, TO = 3) mode of h are shown in Fig. 27 
435M0 (/low 



For M 



lO.OlHz), the match for this is 
7W33 = 0.986921. 



(31) 



As was seen in the g = 1/10 case, we have a somewhat 
large mismatch for the sub-dominant {i — 2, m = 1) and 
(£ = 3, TO = 3) modes. 



D. Amplitude differences 

Here, we focus on the amplitude of the (^ = 2, to = 2) 
mode of h obtained using NR and the SRWZ formalism. 
In Fig. [28] we show the amplitude difference between the 
NR and SRWZ waveforms for the q — 1/10 case. There 
is a large (^ 10%) difference in the amplitudes both in 
the inspiral and merger phases. 

To understand this difference, we compare the ampli- 
tude differences between the NR and SRWZ waveforms 
for the other mass ratio cases. Using the following defi- 
nition for the amplitude difference. 



SA22 = I A 



(SRWZ) 



-A 



(NR), 



(32) 



we find that the amplitude differences for the q = 1/10 
and g 1/15 cases around the maximum amplitude obey 



SA 



(9=1/10) 
22 



2.3 X SA 



(9=1/15) 



1.41 



2.42 



22 
X SA 



(9=1/15) 



(33) 



where 1.41 is the ratio of the symmetric mass ratios 77. 
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FIG. 25: The amplitude (TOP) and phase (BOTTOM) of 
the {£ = 2, m = 2) mode of h for the q = 1/15 case. The 
waveforms were translated in time such that the maximum in 
the amplitude occurs at t = and the phases were adjusted 
by a constant offset such that (j> = at t = — 600M. The 
solid (black) and dashed (red) curves show the NR and SRWZ 
waveforms, respectively. 
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FIG. 26: The amplitude (TOP) and phase (BOTTOM) of 
the {£ = 2, m = 1) mode of h for the q — 1/15 case. The 
waveforms were translated in time such that the maximum in 
the amplitude occurs sd, t — and the phases were adjusted 
by a constant offset such that = at i = — 600M. The 
solid (black) and dashed (red) curves show the NR and SRWZ 
waveforms, respectively. 
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Figure 29 shows the amplitude differences (note that 
t = denotes the time of the maximum amplitude) for 
q ~ 1/10 and q = 1/15. We see that the amplitude 
difference before merger is proportional to the ratio of 
the symmetric mass ratio. This behavior is different from 
our previous analysis in Ref. [H] , where we directly used 
the NR trajectories in the perturbative calculation of the 
waveforms, and found that the amplitude difference in 
the inspiral phase scales like . Since the fitting formula 
is based on the quasicircular evolution, any eccentricity 
of the orbit may create amplitude differences during the 
inspiral phase. 



Figure 30 shows the amplitude differences both for the 
q~ 1/10 and q — 1/100 waveforms. Near the maximum 
amplitude, we find that the amplitude differences scale 
like 



21.5 X <54r'^'™^ 



(34) 



where 8.43 is the ratio of the symmetric mass ratio. 
Note that this rescaling does not capture the behavior 



of the amplitude differences during the inspiral phase. 
Although it is difhcult to obtain any exact relation be- 
tween the amplitude differences during the inspiral (pos- 
sibly due to effects of eccentricity), we do observe non- 
linear effects in the amplitude difference between the NR 
and SRWZ waveforms around the maximum amplitude, 
which we infer from the fact that the amplitude difference 
scales nonlinearly with the mass ratio. 

Figure [3T] shows the amplitude differences both for the 
q = 1/15 and q = 1/100 cases. Near the maximum am- 
plitude, we find that the amplitude differences scale like 



9.5 X M^r'^'™^ 
5.981-26 X M^r'^'™^ 



(35) 



where 5.98 is the ratio of the symmetric mass ratios. 
Again, we see the similar behavior to q = 1/10 and 
q = 1/100 comparison. 

We note that the C transformation in Eq. ( 13 ) also 
affects the amplitude difference. As seen in Fig. 32] if 
we directly use the orbital radius -Rnr to calculate the 
waveforms, we obtain an even larger amplitude than that 
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FIG. 27: The amplitude (TOP) and phase (BOTTOM) of 
the {£ = 3, m = 3) mode of h for the q = 1/15 case. The 
waveforms were translated in time such that the maximum in 
the amplitude occurs at t = and the phases were adjusted 
by a constant offset such that (j> = at t — — 600M. The 
solid (black) and dashed (red) curves show the NR and SRWZ 
waveforms, respectively. 
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FIG. 29: The amplitude difference SA22 of the (^ = 2, m = 2) 
mode of h for the q = 1/10 and 1/15 cases. 
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FIG. 30: The amplitude difference 5A22 of the {I = 2, m = 2) 
mode of h for the q = 1/10 and 1/100 cases. 
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FIG. 28: The amplitude difference of the (^ = 2, m = 2) 
mode of h for the q = 1/10 case. Here, the difference is 
normalized by the NR amplitude. We see ~ 10% difference 
in the inspiral and merger phases. 
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FIG. 31: The amplitude difference &A22 of the (^ = 2, m = 2) 
mode of h for the g = 1/15 and 1/100 cases. 
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obtained using the corrected orbital radius i?Log- This 
behavior is expected based on the fact that while the 
orbital radius is larger, the orbital frequency remains the 



same. Although the orbital radius i?NR gives a similar 
amplitude to the NR waveform in the initial inspiral part, 
we have a much larger amplitude in the merger phase. 
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FIG. 32: The amplitude of the = 2, m = 2) mode of h for 
the q = 1/10 case. We have set the maximum amphtudes at 
i = by the time translation. The solid (black), dashed (red) 
and dotted (blue) curves show the waveforms from the NR, 
SRWZ and SRWZ without the C transformation in Eq. 
respectively. 



and go through the merger of the holes until the forma- 
tion of a common horizon. The fitting formulae Eqs. ([s]) 
and with Eqs. (ITT]) should provide a good approx- 



— NR 

- - SRWZ 

. . SRWZ (without C) 




VI. DISCUSSION 

In the first paper of this series [21] we established that 
perturbation theory of black holes in the extreme-mass- 
ratio expansion q = mi/m2 <C 1, as in Eq. ([6]), can 
faithfully represent an approximation to the full numer- 
ical simulations, if we provide the full numerical relative 
trajectory of the holes as a function of time. In this 
second paper we proceed to model these trajectories us- 
ing full numerical simulations and a parametrized exten- 
sion of 3.5 PN quasi-adiabatic evolutions, Eqs. ([s]) and 



(12 1. This form of the expansion allows us to incorpo- 



rate the PN trajectory at large distances and smoothly 
fit to the full numerical ones up to separations as small 
as i?sch — 3M, at which point we can describe the sub- 
sequent trajectory by a Schwarzschild plunging geodesic. 
In this approximation, we set the values of the energy E 
and angular momentum L of the geodesic to those of the 
previously fitted trajectory evaluated at i?sch = 3M. 

Once we verified the procedure for each of the test 
runs, i.e. for q = 1/10, 1/15 and 1/100, we obtain an 
expression for the fitting parameters as a function of q, as 
given in Eqs. ( 11 ). The other piece of information needed 



by the perturbative approach, in addition to the relative 
trajectory, are the mass and spin of the final black hole 
formed after merger to provide the background metric. 
The spin and mass parameters for the SRWZ equation 
^ can be estimated from empirical formulae ^55, ,56J. 
Note that this formula correctly predicted the remnant 
black hole's mass and spin for the three cases of mass 
ratio studied here. 

The current study can be considered a successful proof 
of principle, and needs to be supplemented by a thor- 
ough set of runs in the intermediate mass ratio regime 
0.01 ^ g ^ 0.1 that start at separations where the post- 
Newtonian approximation is still valid, i.e. R ^ lOAf 



imation within the fitting interval, but extrapolation to 
both the comparable masses or extreme mass ratio could 
quickly lose accuracy. The fitted tracks from our work 
can be compared to alternative approaches to the prob- 
lem such as the EOB [57j and EOBNR [SH HH] and 
the final waveforms to direct phenomenological ones as 
given in Refs. [60l[6T]. Our approach provides a comple- 
mentary model for the intermediate mass- ratio regime 
to the above approaches, which mostly apply to the 
comparable-mass regime. 

The extension of our approach to initially highly- 
spinning black holes can be done essentially through the 
same lines depicted here. Numerical simulations of highly 
spinning holes with small mass ratio are feasible [62] , and 
the trajectories would already include the effect of the 
spins of both holes. The perturbative waveform calcula- 
tion would use the Teukolsky formalism [53] and its sec- 
ond order extension |3I]. Fitting formulae based on the 
extension of spinning post-Newtonian trajectories with 
free parameter dependence on the mass ratio and effec- 
tive spins should be used together with the final matching 
to plunging geodesies on an spinning background. The 
final remnant formulae for spinning binaries have already 
been proposed and fitted to available simulations ^55) . , 56 ] . 
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Appendix A: Perturbative Equations in 1+log 
"trumpet" coordinates 

In this work we used standard Schwarzschild coordi- 
nates to describe the background geometry and trajec- 
tory of the particle. To obtain the latter, we transformed 
the particle's location from the NR coordinates (assum- 
ing they were represented by 1 -I- log trumpet coordi- 
nates of a Schwarzschild spacetime) into the standard 
Schwarzschild coordinates. It is also useful to see the 
form of the perturbative equations in these 1 -t- log trum- 
pet coordinates. This will help to assess potential regions 
where this approach might have limitations. In partic- 
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ular, since the trumpet coordinates are obtained in the 
late-time, stationary limit, wherever there still are strong 
dynamics, the approach could be inaccurate. 

The generalized Regge-Wheeler-Zerilli equations for an 
arbitrary spherically symmetric background were found 
in Ref. |64|- Here, for the sake of simplicity, we use 
only the generalized Regge- Wheeler equation given in 
Ref. '65\. Formally, this equation can be written as 

^RW = Ci + C2 ^RW + C3 ^RW 

+C4 $ V + C5 $w , (Al) 

where the dot and dash denote the time and radial deriva- 
tives, and the coefficients q {i — 1, 2, 3, 4, 5) are deter- 
mined by the background metric. 

For the background metric, we use the notation of [52] . 
which has the form 



~{a^ - (3^)dt'^ + j-dtdr 
+ ^dr^ +R\de^ + sin' 

r 



(A2) 



where a, /3, / and R are functions of the radial coordinate 
r only. The coefficients Ci are then given by 



ci - 2/(r)/?(r), 

C3 = -(/3(r)a'(r)-a(r)/3'(r)) 



a{r) ' 



C4 



a{r) 



(-2 a{r) /3(r) /3'{r) + a{ry a' (r) 



+I3{rf a'{r)) + /(r) /'(r) {a{rf - P{rf) , 
C5 = ^a{TfV{r), (A3) 



where V denotes the Regge- Wheeler potential, 



V{r) = 



1 



R{rY 



6M 
W) 



(A4) 



When using the standard Schwarzschild coordinates, the 
above equation becomes the standard Regge- Wheeler 
equation. We note that for the Zerilli equation, we may 

set m 



V{r) = 



[\^{\ + X)R{rf 



R{rf{XR{r)+?,MY 
+'i\^MR{rf + Q\M^R{r)+mi^] , (A5) 



where A = (£ - 1)(^ -I- 2)/2. 

We now consider the stationary 1 + log slices of the 
Schwarzschild spacetime in isotropic coordinates as an 
approximation to the NR coordinate system. As dis- 
cussed in Ref. [52] . the metric in this coordinate system 
can be obtained numerically. 

Since the Schwarzschild metric is a solution of the vac- 
uum General Relativity field equations, and since we are 
using a Killing lapse and shift [55] , we have 



a(r) = /(r)i?'(r) 



a{rf-li{rf = 1 



2M 

W) 



(A6) 



From the stationary 1 -|- log slice condition, we have 
Ce« = [a" -1+'^-^\r\ (A7) 



where the constant C is determined by requiring that a' 
be regular at the "critical" point (where the denominator 
vanishes) . 



C 

M4 



1 

128' 



(3 + VIO)-" exp(3 - VIO) 



In isotropic coordinates we obtain 

ri?'(r) 



a{r) 



R{r) 



(A8) 



(A9) 



Since the above equations are transcendental, we have 
to find numerical solutions. However, it is also useful to 
study their asymptotic behavior for large r. At large r 
we have, 

R{t) = r + M+- — ---^+0(l/r4), 



a(r) 



1 



M 



4 r 8 r3 

1 NP 1 Af3 

2 ^2 4 

1 C \ AP 

2 M4 



+ 0(l/r^), (AlO) 



and the other two functions j3 and / can be found from 
Eqs. (lAel). 



Inserting these functions into Eqs. (A3), we find that 
the coefficients Ci are given by 
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Cl 



+ 0(l/r3) 



C2 = 
C5 = 



4M 17 M2 

T 7^ 



0(l/r4 




13 A/'' 
,.3 



259 1 C \ , 



39 Af3 



259 3 C \ M4 , 



2 A-/ 17 Af2 

^(^+1) 6M{l+£{e+l)) 1 M2 (84 + 37^(^+1)) 1 Af^ (303 + 79^ + 1)) 



(All) 



r 



If we set C — 0, the above coefRcients are the same 
as those derived in the standard isotropic coordinates 
('^Sch = nso(l + A<f/ (2riso))^) for the Schwarzschild space- 
time. 

An important piece of information in setting up the 
grid refinements, particularly for the smaller black hole 
is the perturbative potential |16| . Using the numerical 
method given in Ref. [52], we can obtain the potential 
term C5. The "Regge- Wheeler" potentials for the stan- 
dard Schwarzschild ("rsch), isotropic (hso), and the trum- 
pet radial coordinates (^Log) (associated to the NR co- 
ordinates) are shown in Fig. 
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while the radial deriva- 
tives of these potentials are shown in Fig. [34] In the 
trumpet coordinates, the radial coordinate covers up to 
''Log = 0, and a and R around rLog = have the following 
forms [52]. 



al^Log) r 



1.091 
Log 



1.092 



(A12) 



where Rq ~ 1.312Af. From the above equation, we see 
that the sphere corresponding to rLog = has a finite 
circumference. 

We observe that the differences between the "trumpet" 
and isotropic coordinates (outside the horizon), as deter- 
mined by differences in the potentials, occurs in a region 
between the horizon and the maximum of the potential. 
This lends credence to the postulate that the main differ- 
ences between the numerical coordinates and the "trum- 
pet" coordinates lie in this region, as well. Indeed, as the 
two BHs approach each other, the conformal function W 
near the larger BH must change (since = at the lo- 
cation of the smaller BH, and varies smoothly outside the 
two punctures). Thus, as the small BH falls through the 
maximum of the potential, the coordinate system must 
become distorted in this important region. Therefore, 
we may expect that the assumption that the trajectory 
in the background spacetime is well described by assum- 
ing the coordinate radius is the trumpet radial coordinate 



breaks down. This may explain the need to match tra- 
jectories to plunging geodesies with slightly larger E and 
L values than one might expect (i.e. E > 1). 

Another use of the perturbative potentials shown in 
Fig. [33 and, in particular, its radial derivative shown 
in Fig. 34J is to guide the set up of the numerical grid 
structure near the black holes [16]. To set up the mesh 
refinement levels, we chose grids that cover the region 
between the zeros of the derivative of the potential. The 
idea is that we need to model the curvature and the 
gravitational radiation emitted by the small BH (which 
drives the merger, and hence the physics). At the ze- 
ros of the derivative of the potentials, the variations are 
minimized. Furthermore the separations between zeros 
increases, naturally leading to a choice of small-width, 
high resolution grids between the first zeros, one step 
lower in resolution between the second two, followed by 
a sequence of coarser grids. 

In order to understand one of the main features of the 
wave propagation, we also calculate the ingoing and out- 
going characteristic speeds for this system [66]. Using 
our definition of the metric, these are given by 



/(r) (a(r)-/3(r)) , 
-fir) (a(r) + /3(r)) 



(A13) 



Figure [35] shows the radial dependence in the trum- 
pet coordinates. Note that inside the horizon, (^Log — 
0.8304Af), both speeds remain negative and that the out- 
going characteristic speed vanishes precisely at the hori- 
zon. Both modes reach their asymptotic speeds (±1), at 
much larger radii. 

It is interesting to observe the differences between the 
usual isotropic coordinates in the standard Schwarzschild 
slicing and the "trumpet" coordinates. The initial data 
are in isotropic coordinates, therefore these differences 
are indicative of how the background has to change from 
the initial slice to the eventual "trumpet" slice. 
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FIG. 33: The solid (black), dashed (blue), and dotted (red) 
curves show the standard Regge-Wheeler (TOP) and Zerilli 
(BOTTOM) potentials in the standard the Schwarzschild co- 
ordinates, isotropic coordinates, and "trumpet" coordinates 
(which are the stationary 1 + log slices of the Schwarzschild 
spacetime in isotropic coordinates), respectively, for the 1 = 2 
mode. Note that the location of the horizon is different in each 
system, rsch < 2M, riao < M/2, rLog < 0.8304A/ are inside 
the horizon. Note that the potential is finite at the horizon, 
and is always positive. 
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